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1 Introduction 

The inverse problem consisting in extracting a source 
map from MEG data is well known to be ill-posed, 
having no unique solution and being numerically in¬ 
stable. Moreover, in the distributed source model, 
first introduced by Hamalainen et al [1], thousands of 
sources are to be reconstructed from only about 200 
simultaneous data samples. 

To solve this problem, a first approach consists in ite¬ 
ratively reducing the search space around emerging 
active areas: Srebro for instance proposed an itera¬ 
tive scheme to limit the head volume to be searched 
at each step [2], We propose here a multiresolution 
iterative procedure in order to deal at each step with a 
constant number of dipoles distributed in a restricted 
area but at an increasing spatial resolution. As an 
extension to Srebro’s work, we build several ellip¬ 
soids instead of a unique area around all the activ¬ 
ity centers and define the new source space at their 
intersection with the cortical surface. We also ex¬ 
tend a more sophisticated regularization strategy than 
the one used by Srebro: our regularization procedure 
extends Baillet et al’s method [3] to multiresolution 
source imaging with enhanced convergence proper¬ 
ties to a unique solution. 

The article is divided as follows: in section 2.1, we 
describe the multiresolution process. In section 2.2 
we present the regularization technique that was used 
and the realistic constraints integrated in the mul¬ 
tiresolution scheme. We finally present some results 
in Section 3 obtained from real MEG data collected 
from a somatosensory experiment with a healthy hu¬ 
man subject. The aim was to reconstruct the under¬ 
lying activity and to compare the performances of the 
multiresolution method to some alternative methods. 

2 Methods 

2,1 Multiresolution procedure 

Let us write the inverse problem of recovering source 
magnitudes from MEG measurements in a linear way: 


M n = GJ n + b n (1) 

where M n is a Nm column vector containing MEG 
measurements at time n, G is the Nm x N gain ma¬ 
trix, J n is an N column vector containing the n th time 
sample of dipole magnitudes and bn is an additive 
noise process. Solving the inverse problem consists 
in obtaining an estimate J n of J n from the data vec¬ 
tor M n . 

Cortical distributed methods aim at determining the 
magnitude of dipoles, perpendicular to the cortical 
surface and located at each node of a mesh repre¬ 
senting the cortical surface. This mesh is obtained 
through a segmentation of structural MRI and gene¬ 
rally contains a large number of nodes (~ 10 4 ). Thus, 
despite their elegant formalism, distributed methods 
are numerically very demanding because of the size 
of J n to be estimated, when such a realistic anatomy 
description is used. In order to overcome this dif¬ 
ficulty, we introduce the following multiresolution 
scheme for iterative distributed source estimation. 

At every resolution but the highest, the number of 
source candidates is limited to N’, a number much 
smaller than the total number of possible sources de¬ 
fined as the nodes of the high-resolution tessellation 
of the cortical surface. The iterative procedure pro¬ 
ceeds to the successive definition of subspaces of lo¬ 
cally increasing source density based on the source 
intensity distribution, estimated at the previous step. 
The iterative procedure runs as follows: 

1. Initialization: source amplitude estimation at the 
lowest resolution 

2. Definition of the new source space: 

(a) Selection of the source with largest ampli¬ 
tude as the centroid of a new search space; 

(b) Estimation of the spatial extension of the 
activity about this centroid as an ellipsoid 
of interest; 

(c) Proceed to the next largest source above a 
threshold left in the source space and repeat 
2.(b); 



3. Definition of the new source space at higher 
resolution within the ellipsoid(s) of interest; 

4. Source amplitude estimation limited to the 
source space defined in 3. using some regula¬ 
rization adapted to the resolution; 

5. Repeat 2.,3.,4. until amplitude estimation is 
done at the highest resolution. 


Let us now describe these different steps. 

Step 2 - At iteration k, step 2 consists in the estima¬ 
tion of both the number of possible active regions and 
their spatial extension around a centroid of maximum 
of intensity. Following the selection of the source 
of largest amplitude, as found at iteration k — 1, the 
extension of the local current density around this lo¬ 
cation is estimated via the definition of an ellipsoid of 
interest which parameters are set following a singu¬ 
lar value decomposition (SVD) of the N x 3 J*~ 1 - 
weighted source coordinates: 

A, = A J.X, (2) 


with 


A J = 


diag(\J% 1 |) 


( 3 ) 


where ,/„ 1 refers to the estimate of J n at iteration 
k — 1 and diag(J^~ x ) is a diagonal matrix which 
non-zero elements are the components of |J^ _1 1. The 
ellipsoid hemi-diameters are set to the obtained singu¬ 
lar values. 

This procedure is repeated starting from the source 
location with the next largest amplitude above a 
threshold //, proportional to the standard deviation 
above the mean value of the distribution J^ _1 in 
the patch, and that does not belong to the previous 
volumes of interest. 

Step 3 - Once the volumes of interest are defined, the 
new source distribution at higher resolution is derived 
from a uniform subsampling of the high resolution 
cortical surface at the intersection with the ellipsoids 
of interest so that the total number of sources reaches 
N'. 


At intermediate resolution, the dipoles can be far 
from each other and represent the average activity 
of a whole region of interest. In this case, regional 
average activity is modeled by three perpendicular 
current dipoles. This holds true until maximal reso¬ 
lution is reached, where there is one dipole per mesh 
node. The normals to the cortical surface can then be 
used as a constraint for the dipole orientation. 


Step 4 - The amplitude estimation is performed with 
a regularization technique described in the following 
section. 

2,2 Estimation procedure 

The regularization technique consists in finding the 
best vector J n according to a given criterion to be 
minimized. This function is a sum of two terms, a 
least-square error term corresponding to the fidelity 
to data and a penalization term L(J n ), corresponding 
to prior knowledge on the source distribution: 

J n = argmin(||M„ - GJ n f + A L(J n )), (4) 

A is a positive scalar that balances the respective con¬ 
tributions to U(J n ) of the data attachment term and 
the prior term L(J n ). A priori information on the 
source distribution to be estimated is introduced in 
order to recover focal source clusters in the source 
image. A first assumption holds that the source 
magnitude pattern is made of areas with smooth inten¬ 
sity changes that may be separated by higher jumps 
in source amplitude: this situation occurs for instance 
between adjacent but functionally non-related corti¬ 
cal areas, e.g. the ones on both sides of a sulcus. 
Moreover, it can be safely considered that relevant 
frequencies of massive neural electrical activity are 
inferior to 100 Hz. Since measurements are oversam¬ 
pled, it is also assumed that dipole magnitude evolves 
smoothly with time. As previously explained, only 
when sources are closer to each other, i.e. at the 
maximal resolution, precise anatomical assumptions 
on the source space are acceptable. Therefore, two 
successive strategies are adopted, one at intermediate 
resolutions and the second one at maximal resolution. 
Let explicit the prior term as: 

L(J n ) = Ls(J n ) + Lt{J n ), (5) 

where L s stands for the spatial prior term and Z* for 
the temporal prior term. 

At each resolution, a temporal neighborhood con¬ 
taining the time sample preceding each source is 
designed. A quadratic cost function is used for the 
temporal prior term. Assuming that J n -\ and J n are 
similar to each other, the orthogonal projection of J n 
on the hyperplane perpendicular to J„_i is "small". 
Thus, L t (J n ) can be written as: 

Z 4 (J„) = ||P„_ 1 J„|| 2 , (6) 

where P n ~\ is the projector onto the hyperplane per¬ 
pendicular to J n -i- 



At intermediate resolution, L s is written as a sum of 
potentials: 

£.(•/*) = 4 . (?) 

where is the i th component of vector J n . 
Non-quadratic (^-functions allow modulated penaliza¬ 
tion of high amplitude sources. Indeed, (j) is growing 
slower than a quadratic function for large amplitudes, 
whereas it has a quadratic behaviour with small am¬ 
plitudes. We chose a convex function to ensure con¬ 
vergence: 

<f>(u) = 2yf 1 + u 2 .K 2 - 2, (8) 

For small values of u, the cost is quadratic whereas 
large values are associated to a linear cost depending 
on the value of K used as a discontinuity threshold. 
Although K could be chosen locally according to 
anatomical priors, K is set constant at intermediate 
resolution for each source, since prior knowledge on 
the activation location is not available. 

At maximal resolution, sources are close to each other 
and relevant anatomy based constraints can be intro¬ 
duced. The idea is to correlate sources that belong to 
a same functional area, and on the opposite to allow 
discontinuity between source amplitudes belonging 
to different functional zones. For this purpose, (j) 
is applied to source gradients. A spatial neighbor¬ 
hood system is designed, with the closest nodes on 
the cortex surface mesh and L s is then written as a 
sum of locally calculated potentials: 

w = EDm), (9) 

i= 1 v=l 

V*" is the gradient operator on sources i and v am¬ 
plitudes and JV„ is the number of neighbors of each 
dipole. In this case K is chosen locally in accor¬ 
dance with anatomical priors, depending on whether 
as K -indexed potential functions <P V are applied to 
gradients between a priori correlated dipoles or not. 
In summary, we adopted two regularization strate¬ 
gies: one at intermediate resolution, and the other at 
maximal resolution where the potential functions are 
applied on gradients to take full advantage of anatom¬ 
ical information. In both cases the (^-function defined 
by equation(8) is used. These non-quadratic criteria 
are minimized by a deterministic descent method, and 
convergence is granted by the convexity of the crite¬ 
ria. Initialization is done with a minimum norm so¬ 
lution, while an empirical adjustment of the hyper¬ 
parameter A is realized. 


3 Results 

In this experiment, we recorded the magnetic brain 
response of an healthy subject to an electric stim¬ 
ulation of the digits (CTF MEG system with 151 
channels, Paris). The purpose was to map the cor¬ 
tical representation of the fingers of right hand of 
the subject and to compare it with prior anatomi¬ 
cal knowledge about somatotopy. We confronted the 
results of the proposed method to those obtained with 
classical approaches : a dipolar method and a min¬ 
imum norm (MN) integrated and not, to a multires¬ 
olution scheme. The results obtained with a min¬ 
imum norm constraint, without multiresolution are 
presented in table 1. The error e mm is the distance be¬ 
tween the source of maximal amplitude and the dipole 
identified by the dipolar method that was used, Dipol- 
eFit software(DF, CTF, Vancouver, Canada). Em is 
the percentage of the source distribution energy con¬ 
tained in the maximal amplitude source. We also give 
the number of sources that had significant amplitude, 
i.e. an amplitude two standard deviations above that 
of the reconstructed distribution. In this case more 
than 9 cm were found between the dipole identified by 
dipolar fit and the maximal amplitude source. Clearly 
these results are unconsistent. But when minimum 
norm regularization was integrated into a multireso¬ 
lution procedure, consistent with DipoleFit’s results 
were obtained, see table 2. The interpretation of the 
obtained distribution is nonetheless problematic given 
the large number of sources with significant ampli¬ 
tude - 79 on average - and the widespread activity. 
Indeed, the maximal amplitude source contained only 
15% of the total energy on average. Results obtained 
with the multiresolution process combined with non¬ 
quadratic regularization to the DipoleFit results are 
given in table 3. This method yielded 45% of the to¬ 
tal energy in the maximal source and 34 significant 
sources on average. 


Table 1: MN without multiresolution compared to DF 


Digits 

1 

2 

3 

4 

5 

e max( cm ) 

9.2 

9.5 

9.1 

9.3 

9.2 

E m (%) 

1.9 

1.0 

1.87 

1.34 

0.86 

active src 

256 

236 

242 

235 

225 




Table 2: MN with multiresolution compared to DF 


Digits 

1 

2 

3 

4 

5 

e max {cm) 

1.5 

5.7 

3.7 

1.3 

1.0 

E m {%) 

15.2 

8.9 

24.1 

9.8 

18.7 

active src 

79 

136 

18 

52 

112 


Table 3: Nonquadratic regularization with multireso¬ 
lution a compared to DF 


Digits 

1 

2 

3 

4 

5 

€max {cm) 

2.3 

1.3 

0.9 

1.2 

0.9 

Em{%) 

91.4 

15.7 

41.8 

51.3 

18.0 

active src 

12 

54 

27 

30 

46 


4 Discussion 

Distributed source models have proven to be an 
interesting direction for extracting realistic images 
of task-related neuronal networks from MEG data, 
but suffer from large source space dimensions and 
bad conditioning of the linear operator. To counter 
this limitation, we have developed a multiresolu¬ 
tion method, combined with a realistic regularization 
based on the anatomy of the subject. This approach 
was validated in an MEG experiment and confronted 
to a dipolar method, and to a distributed source model 
used with a constraint of minimal norm integrated and 
not to a multiresolution scheme. The results showed 
that the multiresolution procedure always yielded a 
solution that was more focal and realistic than with¬ 
out multiresolution. 

Despite their wide use, dipolar methods require some 
prior knowledge on the number of dipoles to fit, i.e. 
the number of activated areas in the brain. The pres¬ 
ence of multiple local minima in the criteria to be 
minimized during the localization procedure of multi¬ 
ple dipoles also induces sensitivity of dipolar methods 
to initialization. In most cases, prior knowledge on 
activity localization is therefore also required. In con¬ 
clusion, dipolar methods are here presented for com¬ 
parison but cannot be considered as a validation. 
Results obtained with a minimum norm constraint 
and no multiresolution were not realistic with a very 
widespread activity all over the cortex. Only when in¬ 
tegrated into a multiresolution scheme, reconstruction 
was slighly improved but was still widespread around 
the activity focus. Therefore, only some anatomy 


based regularization integrated into a multiresolu¬ 
tion procedure allows focal and comparable results in 
terms of localization to dipolar method, without re¬ 
quiring any prior knowledge on the number of activ¬ 
ity centers and on the activity location. 

An interesting clinical implication of the results from 
the multiresolution method is that, contrary to dipo¬ 
lar methods, activation fields are explicitely reco¬ 
vered along the cortical surface. This is consistent 
with the spatial extension of the functional maps 
drawn by electrophysiological studies performed on 
animals, [6]. 

In conclusion, the multiresolution approach appears 
to be a very promising attempt at upgrading MEG to 
the level of a true brain imaging technique. Future de¬ 
velopments should concentrate on the quantification 
of the spatial extension of activity on the cortical sur¬ 
face. 
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